Предварительный анализ данных

library(ggplot2)
library(dplyr)
library(lubridate)
library(vroom)
library(readxl)
library(emoGG)
library(GGally)
library(psych)
library(psych)
library(nortest)
library(qqplotr)
library(reshape)
library(ppcor)
library(car)
library(lm.beta)
library(ellipse)
library(vegan)
library(olsrr)

Считываем данные.

excel_sheets("D:/data/nov/poverty.xls")
## [1] "Sheet1"
dataf <- read_excel("D:/data/nov/poverty.xls", sheet="Sheet1")

Посмотрим на структуру и первые строчки.

head(dataf)
## # A tibble: 6 x 9
##   ...1     BIRTH DEATH INF_DEAH LIFE_M LIFE_F   GNP GROUP  COUNTRY 
##   <chr>    <dbl> <dbl>    <dbl>  <dbl>  <dbl> <dbl> <chr>  <chr>   
## 1 ALBANIA   24.7   5.7     30.8   69.6   75.5   600 EAST_E ALBANIA 
## 2 BULGARIA  12.5  11.9     14.4   68.3   74.7  2250 EAST_E BULGARIA
## 3 CZECHOSL  13.4  11.7     11.3   71.8   77.7  2980 EAST_E CZECHOSL
## 4 FORMER_E  12    12.4      7.6   69.8   75.9    NA EAST_E FORMER_E
## 5 HUNGARY   11.6  13.4     14.8   65.4   73.8  2780 EAST_E HUNGARY 
## 6 POLAND    14.3  10.2     16     67.2   75.7  1690 EAST_E POLAND
dim(dataf)
## [1] 97  9

Это данные по различным показателям уровня жизни для 97 стран. Оставим только один столбец с названием страны, пропуск значения заменен на NA. Группа - это фактор.

datafm <- dataf[ , -9]
names(datafm)[1] <- "COUNTRY"
datafm$GROUP <- as.factor(datafm$GROUP)
head(datafm)
## # A tibble: 6 x 8
##   COUNTRY  BIRTH DEATH INF_DEAH LIFE_M LIFE_F   GNP GROUP 
##   <chr>    <dbl> <dbl>    <dbl>  <dbl>  <dbl> <dbl> <fct> 
## 1 ALBANIA   24.7   5.7     30.8   69.6   75.5   600 EAST_E
## 2 BULGARIA  12.5  11.9     14.4   68.3   74.7  2250 EAST_E
## 3 CZECHOSL  13.4  11.7     11.3   71.8   77.7  2980 EAST_E
## 4 FORMER_E  12    12.4      7.6   69.8   75.9    NA EAST_E
## 5 HUNGARY   11.6  13.4     14.8   65.4   73.8  2780 EAST_E
## 6 POLAND    14.3  10.2     16     67.2   75.7  1690 EAST_E

COUNTRY - название страны

BIRTH - рождаемость на 1000 человек населения

DEAH - смертность на 1000 населения

INF_DEAH - младенческая смертность в возрасте до 1 года на 1000 населения

LIFE_M - продолжительность жизни мужчин

LIFE_F - продолжительность жизни женщин

GNP - валовой национальный продукт на душу населения в долларах США

GROUP - группа страны (в какой части мира находится)

Построим pairs plot.

ggpairs(datafm, columns = 2:8, diag = list(continuous = "barDiag"))

Есть сильно несимметричные распределения и есть неоднородности (видны несколько облаков точек). Категоризующая переменная, объясняющая эту неоднородность - часть мира, в которой расположена страна. Прологарифмируем GNP и раскрасим по группам. Снова строим pairs plot.

datafm$GNP <- log(datafm$GNP)
ggpairs(datafm, columns = 2:8, aes(color = GROUP, alpha = 0.6), diag = list(continuous = "barDiag"))

Зависимости стали более линейными, а распределения более симметричными.

Разделяем на группы.

data1 <- datafm %>% filter(GROUP == "EAST_E")
data2 <- datafm %>% filter(GROUP == "SOU_A")
data3 <- datafm %>% filter(GROUP == "WEST_E_A")
data4 <- datafm %>% filter(GROUP == "ASIA")
data5 <- datafm %>% filter(GROUP == "MID_EAST")
data6 <- datafm %>% filter(GROUP == "AFRICA")

data7<-data

Посмотрим на descriptive statistics. Она показывает основные характеристики распределений по группам. У этой функции следующие значения:

mean - среднее,

sd - стандартное отклонение,

median - медиана,

trimmed - усеченное среднее,

mad - среднее абсолютное отклонение от медианы,

min - минимум,

max - максимум,

range - размах (разница между максимумом и минимумом),

skew - асиммерия,

kurtosis - эксцесс,

se - стандартная ошибка.

Посмотрим на группу AFRICA, в которой 27 стран.

Коэффициент асимметрии для рождаемости равен -0.94. Значит распределение сильно скошено вправо (также это понятно из того, что среднее меньше медианы). Минимум расположен далеко от среднего, так как это значение маловероятно. Оно находится на расстоянии больше, чем 2 сигма (13.43 > 2*5.69). А вероятность быть на расстоянии от мат.ожидания больше двух сигм равна всего 0.046. Эксцесс равен -0.15, значит пик более пологий, а хвосты тяжелее, чем у нормального распределения.

data_fm <- datafm[ , -1]
data_fm <- data_fm[ , -7]
describeBy(data_fm, datafm$GROUP)
## 
##  Descriptive statistics by group 
## group: AFRICA
##          vars  n  mean    sd median trimmed   mad   min    max range  skew
## BIRTH       1 27 44.53  5.69  46.10   45.03  3.56 31.10  52.20  21.1 -0.94
## DEATH       2 27 14.62  4.80  14.00   14.38  5.49  7.30  25.00  17.7  0.48
## INF_DEAH    3 27 99.79 30.58 103.00   99.83 43.00 49.40 154.00 104.6  0.14
## LIFE_M      4 27 50.64  7.10  50.30   50.57  9.19 38.10  64.90  26.8  0.04
## LIFE_F      5 27 54.14  7.03  53.70   54.26  9.49 41.20  66.40  25.2 -0.10
## GNP         6 27  6.20  1.04   6.04    6.17  0.97  4.38   8.58   4.2  0.28
##          kurtosis   se
## BIRTH       -0.15 1.09
## DEATH       -0.85 0.92
## INF_DEAH    -1.34 5.89
## LIFE_M      -1.00 1.37
## LIFE_F      -1.16 1.35
## GNP         -0.65 0.20
## ------------------------------------------------------------ 
## group: ASIA
##          vars  n  mean    sd median trimmed   mad  min    max  range  skew
## BIRTH       1 17 29.62  8.91  30.50   29.97 12.16 11.7  42.20  30.50 -0.25
## DEATH       2 17 10.22  4.66   8.80   10.01  3.85  4.9  18.70  13.80  0.63
## INF_DEAH    3 17 67.72 51.35  64.00   64.24 59.30  6.1 181.60 175.50  0.57
## LIFE_M      4 17 60.49  8.70  62.50   60.87  7.86 41.0  74.30  33.30 -0.60
## LIFE_F      5 17 63.28 10.62  66.10   63.57  9.79 42.0  80.10  38.10 -0.42
## GNP         6 17  6.49  1.35   6.35    6.40  0.72  4.7   9.56   4.86  1.00
##          kurtosis    se
## BIRTH       -1.07  2.16
## DEATH       -1.21  1.13
## INF_DEAH    -0.88 12.45
## LIFE_M      -0.53  2.11
## LIFE_F      -1.04  2.58
## GNP          0.15  0.33
## ------------------------------------------------------------ 
## group: EAST_E
##          vars  n  mean   sd median trimmed  mad  min  max range  skew kurtosis
## BIRTH       1 11 14.76 3.69  13.60   14.01 1.63 11.6 24.7  13.1  1.69     1.91
## DEATH       2 11 10.55 2.08  10.70   10.78 1.78  5.7 13.4   7.7 -0.86     0.14
## INF_DEAH    3 11 17.37 7.06  14.80   16.97 5.19  7.6 30.8  23.2  0.55    -1.05
## LIFE_M      4 11 67.69 2.14  67.20   67.58 2.08 64.6 71.8   7.2  0.36    -1.08
## LIFE_F      5 11 74.99 1.39  74.80   74.98 1.33 72.4 77.7   5.3  0.06    -0.45
## GNP         6  9  7.48 0.48   7.54    7.48 0.27  6.4  8.0   1.6 -1.03     0.10
##            se
## BIRTH    1.11
## DEATH    0.63
## INF_DEAH 2.13
## LIFE_M   0.65
## LIFE_F   0.42
## GNP      0.16
## ------------------------------------------------------------ 
## group: MID_EAST
##          vars  n  mean    sd median trimmed   mad   min   max range  skew
## BIRTH       1 11 33.90  8.63  31.70   33.89 13.20 22.30  45.6 23.30  0.00
## DEATH       2 11  6.75  2.65   7.60    6.73  1.78  2.20  11.5  9.30 -0.12
## INF_DEAH    3 11 47.58 30.77  44.00   45.07 40.03  9.70 108.1 98.40  0.43
## LIFE_M      4 11 64.82  5.01  63.10   64.81  2.08 55.80  73.9 18.10  0.19
## LIFE_F      5 11 67.86  6.06  67.00   68.23  3.26 55.00  77.4 22.40 -0.31
## GNP         6 11  8.52  0.90   8.56    8.52  1.09  7.12   9.9  2.77  0.03
##          kurtosis   se
## BIRTH       -1.81 2.60
## DEATH       -0.96 0.80
## INF_DEAH    -1.08 9.28
## LIFE_M      -0.77 1.51
## LIFE_F      -0.32 1.83
## GNP         -1.39 0.27
## ------------------------------------------------------------ 
## group: SOU_A
##          vars  n  mean    sd median trimmed   mad  min    max range  skew
## BIRTH       1 12 29.18  7.39  28.45   28.55  6.60 18.0  46.60 28.60  0.70
## DEATH       2 12  9.42  5.51   7.65    8.54  1.93  4.4  23.20 18.80  1.49
## INF_DEAH    3 12 51.33 31.70  42.50   48.78 29.43 17.1 111.00 93.90  0.81
## LIFE_M      4 12 62.71  4.93  63.40   63.31  3.78 51.0  68.40 17.40 -0.97
## LIFE_F      5 12 68.53  5.31  68.05   69.19  2.97 55.4  75.10 19.70 -0.90
## GNP         6 12  7.26  0.66   7.35    7.34  0.69  5.8   7.89  2.09 -0.78
##          kurtosis   se
## BIRTH        0.23 2.13
## DEATH        0.85 1.59
## INF_DEAH    -0.71 9.15
## LIFE_M       0.15 1.42
## LIFE_F       0.56 1.53
## GNP         -0.58 0.19
## ------------------------------------------------------------ 
## group: WEST_E_A
##          vars  n  mean   sd median trimmed  mad  min   max range  skew kurtosis
## BIRTH       1 19 12.85 1.94  13.20   12.81 1.93  9.7 16.70  7.00 -0.02    -1.00
## DEATH       2 19  9.43 1.49   9.40    9.45 1.78  6.7 11.90  5.20 -0.08    -1.12
## INF_DEAH    3 19  7.86 1.87   7.50    7.75 0.74  4.5 13.10  8.60  0.95     1.41
## LIFE_M      4 19 71.50 2.66  72.00   71.60 1.93 65.4 75.90 10.50 -0.84    -0.02
## LIFE_F      5 19 78.18 2.30  78.60   78.31 1.93 72.4 81.80  9.40 -0.88     0.24
## GNP         6 18  9.72 0.42   9.82    9.75 0.28  8.7 10.17  1.47 -1.05    -0.09
##            se
## BIRTH    0.45
## DEATH    0.34
## INF_DEAH 0.43
## LIFE_M   0.61
## LIFE_F   0.53
## GNP      0.10

Минимумы и максимумы, какие конкретно страны.

EAST_E

Рождаемость

Венгрия 11.6 Албания 24.7

Смертность

Албания 5.7 Венгрия 13.4

Смертность младенцев

Неизвестная страна (FORMER_E) 7.6 Албания 30.8

Продолжительность жизни мужчин

СССР 64.6 Чехия 71.8

Продолжительность жизни женщин

Румыния 72.4 Чехия 77.7

ВНП

Албания 1.86 Чехия 2.08


SOU_A

Рождаемость

Уругвай 18 Боливия 46.6

Смертность

Венесуэла 4.4 Мексика 23.2

Смертность младенцев

Чили 17.1 Боливия 111

Продолжительность жизни мужчин

Боливия 51 Уругвай 68.4

Продолжительность жизни женщин

Боливия 55.4 Чили 75.1

ВНП

Албания 1.76 Чехия 2.07


WEST_E_A

Рождаемость

Италия 9.7 США 16.7

Смертность

Япония 6.7 Дания 11.9

Смертность младенцев

Япония 4.5 Португалия 13.1

Продолжительность жизни мужчин

Греция 64.4 Япония 75.9

Продолжительность жизни женщин

Португалия 72.4 Япония 81.8

ВНП

Греция 2.16 Финляндия 2.32


ASIA

Рождаемость

Гонконг 11.7 Бангладеш 42.2

Смертность

Гонконг 4.9 Афганистан 18.7

Смертность младенцев

Гонконг 6.1 Афганистан 182

Продолжительность жизни мужчин

Афганистан 41 Гонконг 74.3

Продолжительность жизни женщин

Афганистан 42 Гонконг 80.1

ВНП

Монголия 1.55 Гонконг 2.26


MID_EAST

Рождаемость

Израиль 22.3 Оман 45.6

Смертность

Кувейт 2.2 Иран 11.5

Смертность младенцев

Израиль 9.7 Иран 108

Продолжительность жизни мужчин

Иран 55.8 Израиль 73.9

Продолжительность жизни женщин

Афганистан 55 Гонконг 77.4

ВНП

Иордания 1.96 Неизвестная страна (UNITED_A) 2.29


AFRICA

Рождаемость

Тунис 31.1 Гана 44.4

Смертность

Тунис 7.3 Намибия 12.1

Смертность младенцев

Египет 49.4 Ливия 82

Продолжительность жизни мужчин

Малави 38.1 Нигерия 48.8

Продолжительность жизни женщин

Малави 41.2Нигерия 52.2

ВНП

Мозамбик 1.48 Нигерия 1.77

О виде распределений и о сравнении распределений

Выясним, близки ли распределения к нормальным. Будем анализировать вид распределений признаков по группам, разделяем по частям мира.

Normal probability plot.

df_new <- datafm %>% 
  as.data.frame(datafm) %>% 
  melt(id.vars = c("COUNTRY", "GROUP"))

df1_new <- data1 %>% 
  as.data.frame(data1) %>% 
  melt(id.vars = c("COUNTRY", "GROUP"))
df2_new <- data2 %>% 
  as.data.frame(data2) %>% 
  melt(id.vars = c("COUNTRY", "GROUP"))
df3_new <- data3 %>% 
  as.data.frame(data3) %>% 
  melt(id.vars = c("COUNTRY", "GROUP"))
df4_new <- data4 %>% 
  as.data.frame(data4) %>% 
  melt(id.vars = c("COUNTRY", "GROUP"))
df5_new <- data5 %>% 
  as.data.frame(data5) %>% 
  melt(id.vars = c("COUNTRY", "GROUP"))
df6_new <- data6 %>% 
  as.data.frame(data6) %>% 
  melt(id.vars = c("COUNTRY", "GROUP"))

Отклонения от прямой линии указывают на отклонения от нормальности.

Далее будем сравнивать распределения признаков у SOU_A и MID_EAST, поэтому построим для них normal probability plot.

Видим, что у распределений Рождаемости, Смертности младенцев для группы MID_EAST хвосты более легкие, чем у нормального распределения. (Отличием от тяжелых является направление отклонения от прямой линии для нескольких первых и нескольких последних точек. Для легких хвостов первые несколько точек показывают отклонение от линии над линией, а последние несколько точек показывают отклонение от прямой линии под линией. Для тяжелых хвостов эта схема обратная.)

Normal probability plot для MID_EAST.

ggplot(df5_new, aes(sample = value)) +
        stat_qq_point(size = 2) +
        facet_wrap(~variable, scales = "free") +
        labs(x = "", y = "") +
        stat_qq_line(color = "darkorchid4")+
        ggtitle("Normal probability plot MID_EAST")

Normal probability plot для SOU_A.

ggplot(df2_new, aes(sample = value)) +
        stat_qq_point(size = 2) +
        facet_wrap(~variable, scales = "free") +
        labs(x = "", y = "") +
        stat_qq_line(color = "darkorchid4")+
        ggtitle("Normal probability plot SOU_A")

Распределения похожи на нормальные.

PP-plot.

То есть опять используем графический метод для оценки того, следует ли набор данных заданному распределению (нормальному).

Отклонения от прямой линии указывают на отклонения от нормального распределения.

PP-plot для SOU_A.

ggplot(df2_new, aes(sample = value)) +
        stat_pp_point(size = 2) +
        facet_wrap(~variable, scales = "free") +
        labs(x = "", y = "") +
        stat_pp_line(color = "blue")+
        ggtitle("PP-plot MID_EAST")

PP-plot для MID_EAST.

ggplot(df5_new, aes(sample = value)) +
        stat_pp_point(size = 2) +
        facet_wrap(~variable, scales = "free") +
        labs(x = "", y = "") +
        stat_pp_line(color = "blue")+
        ggtitle("PP-plot MID_EAST")

Распределения похожи на нормальные.

Проверка на нормальность по критериям.

Будем рассматривать отдельно по группам, в какой части мира расположена страна. Проверяем по критериям Лиллиефорса (модификация критерия Колмогорова-Смирнова), Андерсона-Даринга (критерий типа омега-квадрат, придает хвостам больший вес, чем тест KS), критерий хи-квадрат Пирсона для сложной гипотезы нормальности, Шапиро-Уилка (примерно квадрат корреляции между x и y в normal probability ploy).

Не будем проводить тесты для всех признаков всех групп. Рассматривается по 2 группы для признаков Рождаемость, Смертность младенцев, Продолжительность жизни мужчин.

Рождаемость

WEST_E_A

lillie.test(data3$BIRTH)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  data3$BIRTH
## D = 0.097305, p-value = 0.9069
ad.test(data3$BIRTH)
## 
##  Anderson-Darling normality test
## 
## data:  data3$BIRTH
## A = 0.20873, p-value = 0.84
pearson.test(data3$BIRTH)
## 
##  Pearson chi-square normality test
## 
## data:  data3$BIRTH
## P = 3.4737, p-value = 0.4819
shapiro.test(data3$BIRTH)
## 
##  Shapiro-Wilk normality test
## 
## data:  data3$BIRTH
## W = 0.97024, p-value = 0.7813

Гипотеза не отвергается при стандартных уровнях значимости.

SOU_A

lillie.test(data2$BIRTH)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  data2$BIRTH
## D = 0.17612, p-value = 0.3835
ad.test(data2$BIRTH)
## 
##  Anderson-Darling normality test
## 
## data:  data2$BIRTH
## A = 0.40813, p-value = 0.2916
pearson.test(data2$BIRTH)
## 
##  Pearson chi-square normality test
## 
## data:  data2$BIRTH
## P = 8, p-value = 0.04601
shapiro.test(data2$BIRTH)
## 
##  Shapiro-Wilk normality test
## 
## data:  data2$BIRTH
## W = 0.92767, p-value = 0.356

При уровне значимостии 0.05 гипотеза отвергается по критерию хи-квадрат и не отвергается по остальным критериям.

Смертность младенцев

SOU_A

lillie.test(data2$INF_DEAH)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  data2$INF_DEAH
## D = 0.18964, p-value = 0.273
ad.test(data2$INF_DEAH)
## 
##  Anderson-Darling normality test
## 
## data:  data2$INF_DEAH
## A = 0.65002, p-value = 0.06734
pearson.test(data2$INF_DEAH)
## 
##  Pearson chi-square normality test
## 
## data:  data2$INF_DEAH
## P = 4, p-value = 0.2615
shapiro.test(data2$INF_DEAH)
## 
##  Shapiro-Wilk normality test
## 
## data:  data2$INF_DEAH
## W = 0.8588, p-value = 0.04722

При уровне значимостии 0.05 гипотеза отвергается по критерию Шапиро-Уилка и не отвергается по остальным критериям.

MID_EAST

lillie.test(data5$INF_DEAH)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  data5$INF_DEAH
## D = 0.13094, p-value = 0.8625
ad.test(data5$INF_DEAH)
## 
##  Anderson-Darling normality test
## 
## data:  data5$INF_DEAH
## A = 0.28767, p-value = 0.5492
pearson.test(data5$INF_DEAH)
## 
##  Pearson chi-square normality test
## 
## data:  data5$INF_DEAH
## P = 2.6364, p-value = 0.4512
shapiro.test(data5$INF_DEAH)
## 
##  Shapiro-Wilk normality test
## 
## data:  data5$INF_DEAH
## W = 0.93873, p-value = 0.5057

Гипотеза не отвергается при стандартных уровнях значимости.

Продолжительность жизни мужчин

SOU_A

lillie.test(data2$LIFE_M)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  data2$LIFE_M
## D = 0.20085, p-value = 0.1995
ad.test(data2$LIFE_M)
## 
##  Anderson-Darling normality test
## 
## data:  data2$LIFE_M
## A = 0.45735, p-value = 0.217
pearson.test(data2$LIFE_M)
## 
##  Pearson chi-square normality test
## 
## data:  data2$LIFE_M
## P = 1, p-value = 0.8013
shapiro.test(data2$LIFE_M)
## 
##  Shapiro-Wilk normality test
## 
## data:  data2$LIFE_M
## W = 0.89971, p-value = 0.1572

Гипотеза не отвергается при стандартных уровнях значимости.

MID_EAST

lillie.test(data5$LIFE_M)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  data5$LIFE_M
## D = 0.18543, p-value = 0.3617
ad.test(data5$LIFE_M)
## 
##  Anderson-Darling normality test
## 
## data:  data5$LIFE_M
## A = 0.38357, p-value = 0.3309
pearson.test(data5$LIFE_M)
## 
##  Pearson chi-square normality test
## 
## data:  data5$LIFE_M
## P = 2.6364, p-value = 0.4512
shapiro.test(data5$LIFE_M)
## 
##  Shapiro-Wilk normality test
## 
## data:  data5$LIFE_M
## W = 0.94745, p-value = 0.6116

Гипотеза не отвергается при стандартных уровнях значимости.

Как по результатам проверки на нормальность прикинуть, значения асимметрии и эксцесса были около нуля или существенными? Если гипотеза о нормальности отвергается, то значения асимметрии и эксцесса можно считать существенными.

Ящики с усами.

Сравненим распределения в группах с помощью ящиков с усами.

cols <- c("#7edecc", "#79c4db", "#aede81", "#ded17e", "#8d76e8", "#de7eb6")
          

ggplot(datafm, aes(x = GROUP, y = BIRTH, fill = GROUP)) + 
  stat_boxplot(geom = "errorbar",
               width = 0.25) + 
  geom_boxplot(alpha = 0.8,          
               colour = "#474747",
               outlier.colour = 1) + 
  scale_fill_manual(values = cols) +
  ggtitle("Рождаемость")

ggplot(datafm, aes(x = GROUP, y = DEATH, fill = GROUP)) + 
  stat_boxplot(geom = "errorbar",
               width = 0.25) + 
  geom_boxplot(alpha = 0.8,          
               colour = "#474747",
               outlier.colour = 1) + 
  scale_fill_manual(values = cols) +
  ggtitle("Смертность")

ggplot(datafm, aes(x = GROUP, y = INF_DEAH, fill = GROUP)) + 
  stat_boxplot(geom = "errorbar",
               width = 0.25) + 
  geom_boxplot(alpha = 0.8,          
               colour = "#474747",
               outlier.colour = 1) + 
  scale_fill_manual(values = cols) +
  ggtitle("Смертность младенцев")

ggplot(datafm, aes(x = GROUP, y = LIFE_M, fill = GROUP)) + 
  stat_boxplot(geom = "errorbar",
               width = 0.25) + 
  geom_boxplot(alpha = 0.8,          
               colour = "#474747",
               outlier.colour = 1) + 
  scale_fill_manual(values = cols) +
  ggtitle("Продолжительность жизни мужчин")

ggplot(datafm, aes(x = GROUP, y = LIFE_F, fill = GROUP)) + 
  stat_boxplot(geom = "errorbar",
               width = 0.25) + 
  geom_boxplot(alpha = 0.8,          
               colour = "#474747",
               outlier.colour = 1) + 
  scale_fill_manual(values = cols) +
  ggtitle("Продолжительность жизни женщин")

t-критерий.

Смертность младенцев.

Сначала проверим гипотезу о равенстве дисперсий для двух распределений. Рассматриваем Смертность младенцев для SOU_A и MID_EAST.

Критерий Фишера можно использовать только для нормальных распределений, а один из тестов отверг гипотезу о нормальности Смертности младенцев в SOU_A, поэтому будем использовать критерий Левена.

data25 <- datafm %>% filter(GROUP == "SOU_A" | GROUP == "MID_EAST" )
leveneTest(INF_DEAH ~ GROUP, data25, center = "mean")
## Levene's Test for Homogeneity of Variance (center = "mean")
##       Df F value Pr(>F)
## group  1       0  0.998
##       21

Гипотеза не отвергается при стандартных уровнях значимости. Но нельзя считать, что дисперсии равны.

Будем использовать двухвыборочный t-критерий для независимых выборок с равными дисперсиями и с различными дисперсиями, критерий асимптотический. t-критерий точный, когда данные нормальные.

Если бы был сбалансированный дизайн, то не важно равны дисперсии или нет. (Но он не сбалансирован, разные объемы выборки).

t.test(data2$INF_DEAH, data5$INF_DEAH, var.equal=TRUE)
## 
##  Two Sample t-test
## 
## data:  data2$INF_DEAH and data5$INF_DEAH
## t = 0.28688, df = 21, p-value = 0.777
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -23.39176  30.87813
## sample estimates:
## mean of x mean of y 
##  51.32500  47.58182
t.test(data2$INF_DEAH, data5$INF_DEAH, var.equal=FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  data2$INF_DEAH and data5$INF_DEAH
## t = 0.28726, df = 20.921, p-value = 0.7767
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -23.36151  30.84787
## sample estimates:
## mean of x mean of y 
##  51.32500  47.58182

Гипотезы о равенстве математических ожиданий не отвергается при стандартных уровнях значимости.

Продолжительность жизни мужчин.

Аналогично посмотрим на другой признак - Продолжительность жизни мужчин.

Сначала проверяем гипотезу о равенсте дисперсий. Здесь можно использовать критерий Фишера.

var.test(data2$LIFE_M, data5$LIFE_M)
## 
##  F test to compare two variances
## 
## data:  data2$LIFE_M and data5$LIFE_M
## F = 0.9652, num df = 11, denom df = 10, p-value = 0.9474
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2633614 3.4029667
## sample estimates:
## ratio of variances 
##          0.9651967

Гипотеза не отвергается при стандартных уровнях значимости.

Будем использовать двухвыборочный t-критерий для независимых выборок с равными и различными дисперсиями.

t.test(data2$LIFE_M, data5$LIFE_M, var.equal=TRUE)
## 
##  Two Sample t-test
## 
## data:  data2$LIFE_M and data5$LIFE_M
## t = -1.0175, df = 21, p-value = 0.3205
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -6.422153  2.202456
## sample estimates:
## mean of x mean of y 
##  62.70833  64.81818
t.test(data2$LIFE_M, data5$LIFE_M, var.equal=FALSE)
## 
##  Welch Two Sample t-test
## 
## data:  data2$LIFE_M and data5$LIFE_M
## t = -1.0167, df = 20.754, p-value = 0.321
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -6.428751  2.209054
## sample estimates:
## mean of x mean of y 
##  62.70833  64.81818

Гипотезы о равенстве математических ожиданий не отвергается при стандартных уровнях значимости.

Критерий Вилкоксона.

Критерий Вилкоксона непараметрический. Хорош тем, что можно использовать при малых объемах выборки и робастностью (свойство статистического метода, характеризующее независимость влияния выбросов на результат исследования).

Плох тем, что у него мощность меньше, чем у t-критерия.

У нас очень маленькие объемы выборки. Для SOU_A - это 12 индивидов, для MID_EAST - 11 индивидов. Используем критерий Вилкоксона.

wilcox.test(data2$INF_DEAH, data5$INF_DEAH)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data2$INF_DEAH and data5$INF_DEAH
## W = 69.5, p-value = 0.8534
## alternative hypothesis: true location shift is not equal to 0
wilcox.test(data2$LIFE_M, data5$LIFE_M)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  data2$LIFE_M and data5$LIFE_M
## W = 56, p-value = 0.5587
## alternative hypothesis: true location shift is not equal to 0

Гипотезы не отвергаются при стандартных уровнях значимости.

Гипотеза о том, что P(кси1 > кси2) = P(кси1 < кси2). t-критерий, но примененный к рангам.

Критерий Колмогорова Смирнова.

С помощью критерия Колмогорова-Смирнова можно сравнивать распределения в целом (умеет сравнивать формы распределений).

ks.test(data2$INF_DEAH, data5$INF_DEAH)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  data2$INF_DEAH and data5$INF_DEAH
## D = 0.27273, p-value = 0.7868
## alternative hypothesis: two-sided
ks.test(data2$LIFE_M, data5$LIFE_M)
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  data2$LIFE_M and data5$LIFE_M
## D = 0.27273, p-value = 0.7868
## alternative hypothesis: two-sided

Гипотезы о равенстве распределений не отвергаются при стандартных уровнях значимости.

Об анализе зависимостей

У нас неоднородные данные, будем изучать зависимости только внутри групп по-отдельности. Построим сначала pairs plot для группы AFRICA.

Коэффициент корреляции Спирмана не реагирует на монотонные преобразования и почти не реагирует на выбросы. Посмотрим на смертность и продолжительность жизни мужчин, видим, что есть выброс. Значит, коэффициент корреляции Спирмана должен быть больше, чем коэффициент корреляции Пирсона. Проверим это, построив соответствующие матрицы.

ggpairs(data6, columns = 2:8, diag = list(continuous = "barDiag"))

Матрица корреляций Пирсона

dplyr::select(data6, -COUNTRY, -GROUP) %>%
cor(method = "pearson", use = "pairwise.complete.obs") %>%
melt() %>%
ggplot(aes(X1, X2)) +
  geom_raster(aes(fill = value)) +
  geom_text(aes(label = round(value, 3))) +
  scale_fill_gradient2(low=colors()[555], mid=colors()[1], high=colors()[26]) + 
  ggtitle("Pearson AFRICA") +
  theme(axis.text.x = element_text(angle = 50, hjust = 1))

Коэффициент корреляции Пирсона для смертности и продолжительности жизни мужчин равен -0.935.

Матрица корреляций Спирмана

dplyr::select(data6, -COUNTRY, -GROUP) %>%
cor(method = "spearman", use = "pairwise.complete.obs") %>%
melt() %>%
ggplot(aes(X1, X2)) +
  geom_raster(aes(fill = value)) +
  geom_text(aes(label = round(value, 3))) +
  scale_fill_gradient2(low=colors()[555], mid=colors()[1], high=colors()[26]) + 
  ggtitle("Spearman AFRICA") +
  theme(axis.text.x = element_text(angle = 50, hjust = 1))

Коэффициент корреляции Спирмана для смертности и продолжительности жизни мужчин равен -0.951.

Для нормальных распределений коэффициенты Пирсона и Спирмана примерно равны.

Видим достаточно большую корреляцию между продолжительностью жизни (и мужчин, и женщин) и ВНП. Что является причиной, а что следствием? Наверно, чем больше ВНП страны на душу населения, тем больше будет продолжительности жизни людей в этой стране.

Посмотрим на частную корреляцию между Рождаемостью и Смертностью за вычетом ВНП.

(data6 %>%
    dplyr::select(BIRTH, DEATH, GNP) %>%
    pcor(method = "pearson"))$estimate["BIRTH", "DEATH"]
## [1] 0.3536362

Частная корреляция 0.3536 между ними меньше, чем обычная корреляция 0.609.

Посмотрим на частную корреляцию между Продолжительностью жизни мужчин и ВНП за вычетом Смертности.

(data6 %>%
    dplyr::select(LIFE_M, GNP, DEATH) %>%
    pcor(method = "pearson"))$estimate["LIFE_M", "GNP"]
## [1] 0.02495314

Частная корреляция 0.02495 между ними сильно меньше, чем обычная корреляция 0.658.

Линейная регрессия

#data1 <- datafm %>% filter(GROUP == "EAST_E")
data_asia_mid_east <- datafm %>% filter((GROUP == "ASIA")|(GROUP == "MID_EAST"))
data_africa <- datafm %>% filter(GROUP == "AFRICA")

#head(data_asia_mid_east)
ggpairs(data_asia_mid_east, columns = 2:8, diag = list(continuous = "barDiag"))

ggpairs(data_africa, columns = 2:8, diag = list(continuous = "barDiag"))

model_1 <- lm(formula = LIFE_M ~ BIRTH + INF_DEAH + GNP, data = data_africa)
summary(lm.beta(model_1))
## 
## Call:
## lm(formula = LIFE_M ~ BIRTH + INF_DEAH + GNP, data = data_africa)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.1095 -1.9535 -0.2117  1.8078  8.2824 
## 
## Coefficients:
##             Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) 76.04748           NA   10.45968   7.271 2.12e-07 ***
## BIRTH       -0.39253     -0.31450    0.14936  -2.628    0.015 *  
## INF_DEAH    -0.13642     -0.58793    0.02794  -4.883 6.24e-05 ***
## GNP          0.91661      0.13486    0.85326   1.074    0.294    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.398 on 23 degrees of freedom
## Multiple R-squared:  0.7972, Adjusted R-squared:  0.7708 
## F-statistic: 30.14 on 3 and 23 DF,  p-value: 3.821e-08
ellipse_1 <- ellipse(model_1, which=c(2,3), npoints=1000)
#ellipse_1
ggplot() + geom_point(aes(x = ellipse_1[,1], y = ellipse_1[,2]))+xlab("BIRTH")+ylab("INF_DEAH")

ggcorr(data_africa, 
       label = T, 
       label_size = 4,
       label_round = 2,
       hjust = 1,
       size = 5, 
       color = "royalblue",
       layout.exp = 5,
       low = "green3", 
       mid = "gray95", 
       high = "darkorange",
       name = "Correlation")

ols_step_backward_p(model_1)
## [1] "No variables have been removed from the model."
#step(model_1, direction = "backward")
ols_step_forward_p(model_1)
## 
##                             Selection Summary                             
## -------------------------------------------------------------------------
##         Variable                  Adj.                                       
## Step    Entered     R-Square    R-Square     C(p)        AIC        RMSE     
## -------------------------------------------------------------------------
##    1    INF_DEAH      0.6903      0.6779    12.1276    155.7736    4.0274    
##    2    BIRTH         0.7871      0.7693     3.1540    147.6611    3.4085    
##    3    GNP           0.7972      0.7708     4.0000    148.3393    3.3976    
## -------------------------------------------------------------------------
#step(model_1, direction = "forward")
shapiro.test(model_1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model_1$residuals
## W = 0.98271, p-value = 0.9178
ggplot(data_africa ,aes(x = model_1$fitted.values, y = model_1$residuals)) +
  geom_point()+
  xlab("Predicted") +
  ylab("Residuals")

data_africaa <- data_africa[ , -1]
data_africaa <- data_africaa[ , -7]

datt <- data.frame('n' = 1:27, 'country' = data_africa$COUNTRY,
                     'mahalanobis' = mahalanobis(data_africaa, 
                                          colMeans(data_africaa), 
                                          cov(data_africaa)))

ggplot(datt, aes(x = n, y = mahalanobis)) + 
  geom_point(size= 2) + 
  xlab("n") +
  ylab("Mahalanobis' distance")

dat <- data.frame('n' = 1:27, 'country' = data_africa$COUNTRY,
                     'cook' = cooks.distance(model_1))
ggplot(dat, aes(x = n, y = cook)) + 
  geom_point(size= 2) + 
  xlab("n") + 
  ylab("Cook's distance")

ols_plot_resid_stud(model_1)

ols_plot_resid_stud_fit(model_1)

ggplot(data_frame(residuals=rstandard(model_1), studres=studres(model_1)), aes(x=residuals, y=studres))+geom_point()+ geom_abline(scope=1, intercept=0, color = "red")

data_africa_new <- data_africa[ -c(11,15,19),]
model_2 <- lm(formula = LIFE_M ~ BIRTH + INF_DEAH + GNP, data = data_africa_new)
summary(lm.beta(model_2))
## 
## Call:
## lm(formula = LIFE_M ~ BIRTH + INF_DEAH + GNP, data = data_africa_new)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5260 -1.0234  0.1219  1.1299  3.9726 
## 
## Coefficients:
##             Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept) 95.50820           NA    9.48892  10.065 2.84e-09 ***
## BIRTH       -0.49199     -0.37688    0.12434  -3.957 0.000778 ***
## INF_DEAH    -0.17736     -0.77363    0.02377  -7.462 3.36e-07 ***
## GNP         -0.95544     -0.12095    0.85150  -1.122 0.275125    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.578 on 20 degrees of freedom
## Multiple R-squared:  0.885,  Adjusted R-squared:  0.8677 
## F-statistic: 51.29 on 3 and 20 DF,  p-value: 1.421e-09
#new_africa <- data.frame("COUNTRY" = "Land", BIRTH = 38, DEATH = 20, #INF_DEAH = 137, LIFE_M = 60, LIFE_F=70, GNP = 700, GROUP="AFRICA")

new_africa2 <- data.frame(COUNTRY = "Land", BIRTH = 38, INF_DEAH = 137, GNP = 700)
new_africa2
##   COUNTRY BIRTH INF_DEAH GNP
## 1    Land    38      137 700
predict(model_1, newdata = new_africa2, interval="confidence")
##        fit       lwr      upr
## 1 684.0709 -540.6429 1908.785
predict(model_1, newdata = new_africa2, interval="prediction")
##        fit      lwr      upr
## 1 684.0709 -540.663 1908.805